library(readstata13)
library(MASS)
library(foreign)

states <- read.dta13("Barber_Bolton_Thrower_Replication_Data_LSQ.dta")

#first, replicate the main results
summary(model1 <- glm.nb(num_eos ~ as.factor(state_fe) + as.factor(year) + divided_gov1*veto_proof + 
                           house_polarization + rule_review_scale + governor_power + governor_party_pres_vote + 
                           governor_prior_vote + gov_republican + governor_election_year + term_limited + 
                           ln_legislative_staff + ln_legislative_salary + ln_state_income_growth, 
                         data = states[states$gov_republican == 1 | states$gov_democrat == 1,]))
length(model1$resid)

#Run one model leaving out one state at a time - veto-proof majorities
unique <- unique(states$state_fe)
a <- NULL
c <- NULL

for (i in 1:length(unique)){
	
summary(model1 <- glm.nb(num_eos ~ as.factor(state_fe) + as.factor(year) + divided_gov1*veto_proof + 
                             house_polarization + rule_review_scale + governor_power + governor_party_pres_vote + 
                             governor_prior_vote + gov_republican + governor_election_year + term_limited + 
                             ln_legislative_staff + ln_legislative_salary + ln_state_income_growth, 
                           data = states[states$gov_republican == 1 | states$gov_democrat == 1 & states$state %in% unique[-i],]))
  
a <- rbind(a, model1$coef[length(model1$coef)])
b <- sqrt(diag(vcov(model1))[length(model1$coef)])
c <- rbind(c, b)
print(as.character(unique[i]))
}

d <- cbind(a, c)
rownames(d) <- NULL

d <- as.data.frame(d)
names(d) <- c("divided_govt_x_veto_proof", "se")
d$state <- unique

d <- d[order(d$divided_govt_x_veto_proof),]

par(mar = c(6, 4, 4, 2))
plot(1:length(d$divided_govt_x_veto_proof), d$divided_govt_x_veto_proof, ylim = c(-1, 0), pch = 16, xlab = "", ylab = "Coefficient Size", main = "Divided Government, Veto-Proof Majorities, and EO Usage - Omit One State at a Time", axes = F)
segments(x0 = 1:length(d$divided_govt_x_veto_proof), x1 = 1:length(d$divided_govt_x_veto_proof), y0 = d$divided_govt_x_veto_proof - 1.64*d$se, y1 = d$divided_govt_x_veto_proof + 1.64*d$se)
abline(h = 0, lty = 2)

axis(1, at = 1:length(d$divided_govt_x_veto_proof), labels = d$state, las = 2, cex.axis = .7)
axis(2, at = seq(-1, 0, .2), las = 2)
text(20, -.2, "Each Point is Coefficient when Model Run with Labelled State Omitted")
box()

#########################################################

#first, replicate the main results
summary(model1 <- glm.nb(num_eos ~ as.factor(state_fe) + as.factor(year) + divided_gov1*house_polarization + 
                           veto_proof + rule_review_scale + governor_power + governor_party_pres_vote + 
                           governor_prior_vote + gov_republican + governor_election_year + term_limited + 
                           ln_legislative_staff + ln_legislative_salary + ln_state_income_growth, 
                         data = states[states$gov_republican == 1 | states$gov_democrat == 1,]))
length(model1$resid)

#Run one model leaving out one state at a time - polarization
unique <- unique(states$state_fe)
a <- NULL
c <- NULL

for (i in 1:length(unique)){
	
  summary(model1 <- glm.nb(num_eos ~ as.factor(state_fe) + as.factor(year) + divided_gov1*house_polarization + 
                             veto_proof + rule_review_scale + governor_power + governor_party_pres_vote + 
                             governor_prior_vote + gov_republican + governor_election_year + term_limited + 
                             ln_legislative_staff + ln_legislative_salary + ln_state_income_growth, 
                           data = states[states$gov_republican == 1 | states$gov_democrat == 1 & states$state %in% unique[-i],]))
  
a <- rbind(a, model1$coef[length(model1$coef)])

b <- sqrt(diag(vcov(model1))[length(model1$coef)])
c <- rbind(c, b)
print(as.character(unique[i]))
}

d <- cbind(a, c)
rownames(d) <- NULL

d <- as.data.frame(d)
names(d) <- c("divided_govt_x_polariziation", "se")
d$state <- unique

d <- d[order(d$divided_govt_x_polariziation),]

dev.new(height = 5, width = 9)
par(mar = c(6, 4, 4, 2))
plot(1:length(d$divided_govt_x_polariziation), d$divided_govt_x_polariziation, ylim = c(0, .8), pch = 16, xlab = "", ylab = "Coefficient Size", main = "Divided Government, Polarization, and EO Usage - Omit One State at a Time", axes = F)
segments(x0 = 1:length(d$divided_govt_x_polariziation), x1 = 1:length(d$divided_govt_x_polariziation), y0 = d$divided_govt_x_polariziation - 1.64*d$se, y1 = d$divided_govt_x_polariziation + 1.64*d$se)
abline(h = 0, lty = 2)

axis(1, at = 1:length(d$divided_govt_x_polariziation), labels = d$state, las = 2, cex.axis = .7)
axis(2, at = seq(0, 1, .2), las = 2)
text(20, .6, "Each Point is Coefficient when Model Run with Labelled State Omitted")
box()

#########################################################

#first, replicate the main results
summary(model1 <- glm.nb(num_eos ~ as.factor(state_fe) + as.factor(year) + divided_gov1*rule_review_scale + 
                           veto_proof + house_polarization + governor_power + governor_party_pres_vote + 
                           governor_prior_vote + gov_republican + governor_election_year + term_limited + 
                           ln_legislative_staff + ln_legislative_salary + ln_state_income_growth, 
                         data = states[states$gov_republican == 1 | states$gov_democrat == 1,]))
length(model1$resid)

#Run one model leaving out one state at a time - rule review
unique <- unique(states$state_fe)
a <- NULL
c <- NULL

for (i in 1:length(unique)){
  
summary(model1 <- glm.nb(num_eos ~ as.factor(state_fe) + as.factor(year) + divided_gov1*rule_review_scale + 
                        veto_proof + house_polarization + governor_power + governor_party_pres_vote + 
                        governor_prior_vote + gov_republican + governor_election_year + term_limited + 
                        ln_legislative_staff + ln_legislative_salary + ln_state_income_growth, 
                        data = states[states$gov_republican == 1 | states$gov_democrat == 1 & states$state %in% unique[-i],]))
  
a <- rbind(a, model1$coef[length(model1$coef)])
b <- sqrt(diag(vcov(model1))[length(model1$coef)])
c <- rbind(c, b)
print(as.character(unique[i]))
}

d <- cbind(a, c)
rownames(d) <- NULL

d <- as.data.frame(d)
names(d) <- c("divided_govt_x_scope_review", "se")
d$state <- unique

d <- d[order(d$divided_govt_x_scope_review),]

dev.new(height = 5, width = 9)
par(mar = c(6, 4, 4, 2))
plot(1:length(d$divided_govt_x_scope_review), d$divided_govt_x_scope_review, ylim = c(-1, 0), pch = 16, xlab = "", ylab = "Coefficient Size", main = "Divided Government, Rule Review, and EO Usage - Omit One State at a Time", axes = F)
segments(x0 = 1:length(d$divided_govt_x_scope_review), x1 = 1:length(d$divided_govt_x_scope_review), y0 = d$divided_govt_x_scope_review - 1.64*d$se, y1 = d$divided_govt_x_scope_review + 1.64*d$se)
abline(h = 0, lty = 2)

axis(1, at = 1:length(d$divided_govt_x_scope_review), labels = d$state, las = 2, cex.axis = .7)
axis(2, at = seq(-1, 0, .2), las = 2)
text(20, -.8, "Each Point is Coefficient when Model Run with Labelled State Omitted")
box()



